home *** CD-ROM | disk | FTP | other *** search
/ Mac Magazin/MacEasy 21 / Mac Magazin and MacEasy Magazine CD - Issue 21.iso / Wissenschaft & Technik / yorick12vr1-nofpu folder / include / test3.i < prev    next >
Text File  |  1995-07-26  |  4KB  |  142 lines
  1. /*
  2.    TEST3.I
  3.    Pure scalar math problem for timing Yorick.
  4.  
  5.    $Id: test3.i,v 1.1 1993/08/27 18:50:06 munro Exp $
  6.  */
  7. /*    Copyright (c) 1994.  The Regents of the University of California.
  8.                     All rights reserved.  */
  9.  
  10. func test3(npass)
  11. /* DOCUMENT test3
  12.          or test3, npass
  13.      Computes the ratio r which solves 1 + r^2 + r^3 +...+ r^n = s,
  14.      given n and s.  If NPASS is given, the calculation is repeated
  15.      that many times (actually the equation is solved many times for
  16.      each pass).  The worker routine invgeom can actually be
  17.      vectorized; the vector version is gseries_r in series.i.  */
  18. {
  19.   if (is_void(npass) || npass<=0) npass= 1;
  20.  
  21.   extern s, n, answer;
  22.   nsums= numberof(answer);
  23.   r= 0.0*answer;
  24.  
  25.   now= split= array(0.0, 3);
  26.   timer, now;
  27.   m= npass;
  28.   while (m--) {
  29.     for (i=1 ; i<=nsums ; i++) r(i)= invgeom(s(i), n);
  30.   }
  31.   timer, now, split;
  32.   timer_print, "Time per pass", split/npass, "Total time", split;
  33.  
  34.   if (anyof(abs(r-answer)>1.e-9*answer))
  35.     write, "***WARNING*** values returned by invgeom are not correct";
  36. }
  37.  
  38. #if 0
  39. CPU seconds per pass:
  40.             IDL    Yorick     Basis       DAP
  41. HP730        -      0.11       5.08      5.00
  42. Solbourne  0.31     0.21       11.8      13.7
  43.     (varies by ~10%)
  44.  
  45.      forum[10] yorick
  46.       Yorick ready.  For help type 'help'
  47.      > #include "/home/miggle/munro/Yorick/include/test3.i"
  48.      > test3,1
  49.             Timing Category     CPU sec  System sec    Wall sec
  50.               Time per pass       0.190       0.000       0.230
  51.              Total time       0.190       0.000       0.230
  52.      > test3,51
  53.             Timing Category     CPU sec  System sec    Wall sec
  54.               Time per pass       0.206       0.000       0.207
  55.              Total time      10.520       0.000      10.560
  56.  
  57.      forum[6] time idl
  58.      IDL> .rnew /home/miggle/munro/Yorick/include/test3
  59.      % Compiled module: INVGEOM.
  60.      % Compiled module: TEST3.
  61.      % Compiled module: SPAN.
  62.      % Compiled module: $MAIN$.
  63.      IDL> test3,1
  64.      IDL> exit
  65.      0.5u 0.3s 0:08 11% 0+928k 0+2io 0pf+0w
  66.      forum[7] time idl
  67.      IDL> .rnew /home/miggle/munro/Yorick/include/test3
  68.      % Compiled module: INVGEOM.
  69.      % Compiled module: TEST3.
  70.      % Compiled module: SPAN.
  71.      % Compiled module: $MAIN$.
  72.      IDL> test3,51
  73.      IDL> exit
  74.      16.0u 0.4s 0:41 39% 0+1400k 0+0io 0pf+0w
  75.  
  76.      forum[13] basis
  77.      Initializing Basis System
  78.      Basis 7.0
  79.      Initializing EZCURVE/NCAR Graphics
  80.      Ezcurve/NCAR 2.0
  81.      Basis> echo=0
  82.      Basis> read /home/miggle/munro/Yorick/include/test3.bas
  83.      End of input file /home/miggle/munro/Yorick/include/test3.bas
  84.      Resuming input from TERMINAL
  85.      Basis> test3(1)
  86.  
  87.     CPU (sec)   SYS (sec)
  88.        11.817       0.017
  89.  
  90.      forum[14] dap
  91.      DAP> echo=0
  92.      DAP> read /home/miggle/munro/Yorick/include/test3.bas
  93.      DAP> test3(1)
  94.  
  95.     CPU (sec)   SYS (sec)
  96.        13.700       0.033
  97. #endif
  98.  
  99. /* invgeom(s, n)
  100.      returns the ratio r of the finite geometric series, given the sum s:
  101.         1 + r + r^2 + r^3 + ... + r^n = s     */
  102. func invgeom(s, n)
  103. {
  104.   nn= n+1;
  105.   if (nn<2) return s-1.0;
  106.  
  107.   /* compute an approximate result which has exact values and
  108.      derivatives at s==1, s==n, and s->infinity */
  109.  
  110.   /* different approximations apply for s>nn and s<nn */
  111.   if (s>nn) {
  112.     pow= 1.0/n;
  113.     npow= nn^pow - 1.0;
  114.     n2r= 1.0/(nn-2.0);
  115.     A= (2.0-nn*npow)*n2r;
  116.     B= (2.0*npow-nn*pow)*nn*n2r;
  117.     r= s^pow - pow + A*(nn/s)^pow + B/s;
  118.   } else {
  119.     sn= (s-1.0)/n;
  120.     n2r= 1.0/(nn*nn);
  121.     r= 1.0 - 1.0/s + n2r*sn*sn*(nn+1.0 - sn);
  122.   }
  123.  
  124.   /* Polish the approximation using Newton-Raphson iterations.  */
  125.   for (;;) {
  126.     rr= r-1.0;
  127.     rn= r^(nn-1);
  128.     delta= rr*s - (r*rn-1.0);
  129.     if (abs(delta)<=1.e-9*abs(rr*s)) break;
  130.     r+= delta/(nn*rn-s);
  131.   }
  132.   /* try to get it to machine precision */
  133.   if (delta) r+= delta/(nn*rn-s);
  134.  
  135.   return r;
  136. }
  137.  
  138. /* Set up a bunch of values to compute.  */
  139. answer= span(0.0, 2.0, 200);  /* number must be even to avoid 1.0000 */
  140. n= 100;
  141. s= (answer^(n+1)-1.0)/(answer-1.0);
  142.